Bubble generation in a twisted and bent DNA-like model 
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The DNA molecule is modeled by a parabola embedded chain with long-range interactions between 
twisted base pair dipoles. A mechanism for bubble generation is presented and investigated in two 
different configurations. Using random normally distributed initial conditions to simulate thermal 
fluctuations, a relationship between bubble generation, twist and curvature is established. An 
analytical approach supports the numerical results. 

PACS numbers: 05.45. Yv, 63.20.Ry, 63.20.Pw, 87.15.Aa 



O 

OO 

O 

o 



> 
X 



I. INTRODUCTION 

In recent years, biomolecular modeling has received an 
ever increasing amount of attention, especially focused 
on the DNA molecule as well as protein structures. The 
basic structure of DNA is fairly well understood since the 
discovery of Crick and Watson , but it is becoming in- 
creasingly apparent that structure alone does not explain 
its complex functionality sufficiently 0, 0, 13 • 

An example is the mechanism leading to bubble gen- 
eration in DNA, in which the two polypeptide strands 
open to allow replication of the molecule, processing of 
proteins or complete strand separation (denaturation). 
Thermal fluctuations at physiological temperatures and 
nonlinear localizations are expected to produce bubbles 
when geometrical features, such as twist and curvature, 
are taken into account. 

In initial works investigating the denaturation bubble, 
the geometrical features of the molecule were essentially 
neglected and energy localization was mostly attributed 
to inhomogeneities and impurities in the lattice chain, 
which may model the action of an enzyme JE 0j8y9, 
EmiH, or nonlinear excitations H [H Q Cllll 
Il7| . Also, discreteness plays an important role for the 
localization of these excitations. The inhomogeneities 
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have been modeled by different masses at various chain 
sites 0, EH E3 j by conformational defects [T^ l or b y 
changes in the coupling between molecular sites [3, El • 
Also, different on-site potentials |8l ll8lll9ll2CI| have been 
used as inhomogeneities, corresponding to the different 
number of hydrogen bonds between the two strands. In 
DNA, the AT base pair connects through two hydrogen 
bonds, whereas the CG base pair has three. It has to 
some extent been experimentally verified that bubbles 
form at AT-rich sections of the DNA molecule [2l|, but 
recent work 22] also suggests that other mechanisms are 
involved. 

Recently, both long-range dipole-dipole interaction 
[Ulli, helicity [HES and curvature jHIHIHl have 
been included in the nonlinear transport theory, as well 
as combinations of these effects [3(j, bill |33. HE uM ■ ^ has 
been shown that chain geometry induces effects similar 
to those of impurities [271 l28l 123 . l32| . 

In biological environments, thermal fluctuations are al- 
ways present and have been considered in Refs. 0, 
l36l l38l l3^ |. for example. In these references it was 
shown that solitons or discrete breathers can be gener- 
ated from or exist among random thermal fluctuations. 

The aim of the present work is to study an augmented 
Peyrard-Bishop model of the DNA molecule |4jj • We in- 
clude both dipole-dipole long-range interaction and chain 
geometry in the form of a rigid, parabola embedded 
chain. Elaborating on earlier work |41|, we show how 
both chain curvature and twist can initiate bubble gen- 
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FIG. 1: DNA chain embedded on a parabola in the xy plane. 
Sites (i n , y n ) with y„ = (fc/2)(a;„— j) 2 , indicated by light grey 
dots. Base pair dipoles, orthogonal to the parabola, shown as 
dark grey arrows. Curvature n — 2 and twist r = 1. The z 
axis (not shown) forms a right-handed system with x and y 
axes. Left, the on-site case, 7 = 0. Right, the intersite case, 
7 = 1/2. 



eration in the DNA molecule. We consider two differ- 
ent chain configurations and randomly distributed initial 
conditions, modeling physiological temperatures. 

In Sec. [nj we introduce the model Hamiltonian and 
equations of motion and discuss relevant parameter val- 
ues as well as the dipole interaction and chain geometry. 
In Sec. IIIII numerical investigations are performed and 
the results are supported by an analytical approach in 
Sec. IIVI Finally, Sec. contains a summary and a dis- 
cussion. 



previously been used to describe thermal denaturation in 
DNA, see Ref. M. 
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FIG. 2: The Morse potential VWse(un) = (e~"" - l) 2 (solid 
curve). Effective potential, Eq. ©, in the analytical approx- 
imation (dashed curve). 

The LRI is a dipole-dipole interaction, with coefficients 
J nm given by ilE2 



J 71 



J \d>n ' d m 3 {d n • T nra S j (dm • fftm)} 



(2) 



where r n and d n are the position vector and the unit 
dipole vector at the nth site, respectively, and r nrn de- 
notes the unit vector from the nth to the mth site 



II. MODEL 

We consider parabola embedded chains, as illustrated 
in Fig. n] The base pair sites are embedded along a 
parabola in the xy plane at uniform distances. The dipole 
moments of the base pairs, represented by arrows, are or- 
thogonal to the parabola and twisted in the orthogonal 
plane. Within the nth base pair, the deviation from the 
equilibrium transverse distance between the bases is de- 
noted u n |40j. The intrasite dynamics of the base pairs 
is governed by a Morse potential, see Fig. |3 and we as- 
sume a harmonic intersite coupling between neighboring 
base pairs. The coupling — or stacking — parameter, C, 
remains constant along the chain, i.e., independent of 
curvature and twisting. 

As a result, using the scalings and parameter values 
presented below, we obtain the dimensionless Hamilto- 
nian, 



n=N , n 
n=-N v 



U n+1 - U n ) 



+ (e u "-l) 2 + ^^ Jn m u n u m ^, (1) 



where the prime indicates m ^ n in the last summa- 
tion, accounting for the long-range interaction (LRI) be- 
tween the dipoles. The total number of sites are thus 
N T = 2N + 1. Without the LRI, this Hamiltonian has 



(3) 



We note that the geometry of the chain only comes into 
play through the LRFs 0] • 

Dimensionless variables have been introduced as 



t = t/t 



r n /l, and H = H/D, 



where original physical variables are indicated by tildes. 
I is the constant intersite distance between neighboring 
base pairs, a and D are parameters in the dimensional 
Morse potential, D ( e ~ aUn — l) , and the time constant 

is given by to = \/ M / Da 2 , where M is the mass of a 
base pair. The dimensionless parameters C and J are 
given by 



C = C/(Da 2 ) and J = 2J /l 3 Da 2 , 



(4) 



with J = q 2 /47T£ , where q denotes the dipole charge of 
the base pair and £q is the dielectric constant. 



A. Parameter values 

We use the parameter values D — 0.04 eV (= 0.64 x 

10~ 20 J), a = 4.45 A" 1 (= 4.45 x 10 10 m" 1 ), M = 
300 a.m.u. (= 5.00 x 10~ 25 kg) and the coupling param- 
eter C = 0.06 eV/A (= 0.96 J/m 2 ), which have been 
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Symbol 



Parameter 



Physical value 



C 
D 
a 
M 
9 

Jo 
I 

to 



stacking 
Morse depth 



0.96 J/m 2 
0.64 x 10" 20 J 

base pair mass 5.00 x 10 - 25 kg 
dipole charge 1.17 x 10~ 19 C 
interaction strength 0.90 x 10~ 28 Jm 
lattice constant 3.4 x 10 -10 m 
time constant 0.20 x 10~ 12 s 



TABLE I: Physical parameters for the DNA molecule. 

widely used in DNA-like models [HH |13 . The dimen- 
sionless stacking parameter then becomes C — 0.075, 
which we use throughout the paper. We note that C- 
values between 0.003 eV/A 2 and 31.7 eV/A 2 
have been reported in the literature. 

The resulting time constant, to = 0.20 ps, is in the 
picosecond range, as seen in Table [|] 

The dipole moment for the base pair, d^p, is the geo- 
metric sum of the moments for the bases, which range 
between 3 and 7 debye |44|- In our simulations we 
use ddip ~ 7 debye. The corresponding dipole charge 
q = rfdip/adip, where adi P is the equilibrium distance 
between the dipole charges 2 A 23]) then becomes 
1.17 x 10~ 19 C, yielding J = 0.5 0. 

B. Geometry 

We introduce the twist angle, (f> n , which the dipoles 
create with the z direction, in the form of a kink J5|, 
where 7 gives the position of the kink center, 



2 arctan 



-T (n— 7) 



(5) 



Thus, 7 is in the on-site case, but 7 = 1/2 in the 
intersite case. The most important difference between 
the two is that the on-site case always has a dipole aligned 
in the bending plane at n — (</>o = 7r/2) . 

As illustrated in Fig. [3| the larger twist, r, the faster 
the twist angle changes — especially in the region of max- 
imal twist and curvature. Note that the maximal slope 
of the twist function occurs at 7. In the limit r — > 00, 
Eq. JSJ) gives <j> n — ir for n < 0, and <j> = for n > 1 for 
both cases. In the same limit, at n = 0, the on-site case 
has cj>o — 7r/2, whereas the inter-site case has <po = tt. In 
the limit r — > 0, all tj> n = ir/2 in both cases. 

On the parabola embedded chain, y n = 
(k/2) (x n — 7) , the resulting unit dipole vectors 
then become 

d» = (-€ n K(x n - 7)sin0„,£„sin</>„,cos0„), (6) 



with £„ = l/yl + k 2 {x n — 7) 2 (corresponding to the 
dark grey arrows in Fig.^). In the on-site case, :co~7 = 0, 



and the other site positions are numerically calculated to 
fulfill the requirement that the distance between adjacent 
sites is always unity. In the intersite case, we define xq — 
7 = —1/2 and x\ ~ 7 = 1/2 and compute the other 
sites in a similar way. Thus, the axis of symmetry passes 
through the site n = in the on-site, 7 = 0, case and 
passes through the middle of the bond connecting sites 
n = and n = 1 in the intersite, 7 = 1/2, case |47|. See 
Fig.IU 
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FIG. 3: (Color online) The twist angle, < 
of the twist, r, around the chain center, 
intersite case, 7 = 1/2. 



in, for various values 
On-site case, 7 = 0, 



It is difficult to determine the actual value of the twist, 
t. The standard value for the twist angle between neigh- 
boring base pairs in an undisturbed DNA molecule is 
about 36° Hi El 113. Requiring that a twist of 180° 
should be achieved over five sites, thus corresponds to a 
value of t about 1 (Fig.|3) in equilibrium and larger when 
twisted more. 

Similarly, the bending of a DNA chain can be ap- 
proximately determined from experimental results. In 
Refs. |27l l3(il ] a parabolic approximation with curvature 
parameter k < 4 was used. Here, we shall only consider 
the range k < 2 as larger curvature has never been ob- 
served in measurements. 



III. SIMULATION RESULTS 

From the Hamiltonian Q we obtain the equations of 
motion 



u n + C (2u n - U n -l - u n+ i) 

-2e (e l) -\- ^ ^ Jnm^m 



0. 



(7) 



In the following we solve these equations numerically 
using a fourth order Runge-Kutta solver with free bound- 
ary conditions w_tv-i = u_jv and u^v+i = u«- In all 
simulations the relative change of the Hamiltonian is less 
than 10~ 5 and we consider a chain with Nt = 99 sites. 

In our previous work considering an approximate 
dipole-dipole interaction on a wedge shaped chain |4l| 
we investigated random initial conditions. This was done 
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as nonlinear excitations are known to be generated from 
randomness 0, |3^, [33, |23, H|J ■ As the outcome of colli- 
sions between nonlinear excitations depends strongly on 
their relative phases, we made 500 different random ini- 
tial conditions to be able to find effects independent of 
the random phases. The random initial conditions cre- 
ated nonlinear excitations, which led to bubble genera- 
tion at various collision sites. We found that bending of 
the chain caused the bubble generation to localize at the 
bent region. 

In the following we use a similar approach to investi- 
gate the relationship between twist and curvature in this 
realistic expression for the dipole interaction © . We con- 
sider random initial conditions: Initial displacements set 
to zero, i.e., u n (0) = for all n, while initial velocities of 
the chain sites are normally distributed with mean value 
(w„(0)) = and standard deviation Gu n . The standard 
deviation is chosen to be <Ju n = 1.156, corresponding a 
temperature of T w 310 K. 

The system dynamics is simulated for 100 different re- 
alizations of the initial conditions with stepsize 0.01 in 
time. The simulations run for 100 time units (corre- 
sponding to 20 ps) or until the Hamiltonian is no longer 
conserved. As mentioned, the constant dipole-dipole in- 
teraction coefficient J = 0.5 is used in all the simulations. 

For different values of twist or curvature, the same ini- 
tial condition can result in very different behavior. Con- 
sider Fig. ^ where the same initial condition is simulated 
for two different values of the twist, t, in the on-site case. 
With the larger twist (dashed curve) the amplitude grows 
in an exponential-like manner (see Sec. HVfl . 

80 r 

°K 60 
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FIG. 4: Evolution of the center site amplitude for n — 1.0 and 
r = 6 (dashed curve) and r = 4 (solid curve) for an on-site 
case. 



A. Results for the on-site case 

The combined effect of both curvature and twist for 
the on-site case is depicted in Fig. |S| where the shaded 
region corresponds to bubble generation. We see that 
bubbles are generated for strong twist and strong curva- 
ture, which is expected here. For the on-site case, the 
strongest attraction between dipoles occur at the sites 
n = — 1 and n = 1 (which are almost antiparallel for 
strong twist). Increasing the curvature bring these next- 
to-center sites closer, which augments the dipole-dipole 
interaction @. As the inset (a) shows, the amplitude 
increase is localized to the region of maximal twist and 
curvature. We are aware that in reality, whole regions of 
DNA base pairs move apart during denaturation. There- 
fore, our mechanism should only be perceived as a pre- 
cursor for bubble generation. 

The exact shape of Fig. [5] depends on the chosen am- 
plitude threshold — as well as system parameters — but its 
qualitative shape is unchanged. Close to the region bor- 
der, only a few of the simultations results in bubble gen- 
eration, but this number increases as one proceeds in 
the direction of stronger twist and larger curvature (i.e., 
towards the upper right corner). We note that for the 
on-site case the amplitude is increased at the three cen- 
ter sites n = — 1, n = 0, and n = 1 as indicated in the 
inset (a). 




2- 



0.5 ^ 1 1.5 2 

Curvature, k 

FIG. 5: Region of bubble generation in the on-site case. Light 
grey region, bubble generation. Dark grey region, transition 
region indicating the uncertainty of the simulations. White 
region, no bubble generation. Insets show the displacement, 
u n , versus site, n, at the points (a) and (b) at simulation times 
(a) t — 51 and (b) t = 100 for the same initial condition. 



At the end of each simulation, we examine the am- 
plitudes. We use a threshold value of u n = 100, corre- 
sponding to about 20 A, i.e., twice the equilibrium dis- 
tance between base pairs. If the threshold is exceeded in 
at least two adjacent sites, we consider this as a precur- 
sor for bubble generation. In the following sections, we 
have depicted the regions of curvature, n, and twist, r, 
where at least one of the 100 simulations result in bubble 
generation. 



B. Results for the intersite case 

In the intersite case, Fig. the picture is different. 
First of all, the twist needed for bubble generation is 
smaller. Second, the dependence on the curvature is less 
pronounced. This is because the dipoles, that for strong 
twist are antiparallel, in this case are neighbors. Since, 
in the framework of our model, the chain has constant 
distance between adjacent sites, increasing the curvature 
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does not increase the tendency for bubbles to be gener- 
ated. 

In fact, the opposite is the case: As the dipole twist 
is perpendicular to the chain, increasing the curvature 
has the consequence that the center dipoles interact in a 
less attracting way, since the center dipoles become more 
antiparallel (Fig. |SJ. 

Note that in this case, the displacement at the two 
center cites, n = and n — 1, is increased as the inset 
(a) shows. Close to the region border, the number of 
simulations resulting in bubble generation is small, but 
it increases as one increases the twist or decreases the 
curvature (i.e., moves to the upper left corner). 



2.5i — r 
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FIG. 6: Region of bubble generation in the intersite case. 
Light grey region, bubble generation. Dark grey region, tran- 
sition region indicating the uncertainty of the simulations. 
White region, no bubble generation. Insets show the dis- 
placement, u n , versus site, n, at the points (a) and (b) at 
simulation time (a) t = 45 and (b) t = 100 for the same 
initial condition. 



the bending point, there exist an effective potential well 
for t larger than about 0.5. This corresponds to the be- 
havior found in Figs.ElandEl The depth of the potential 
well increases with increasing twist, until a saturation is 
reached at r w 5. The existence of a potential well is a 
necessary, but not sufficient, condition for bubble gener- 
ation. 




FIG. 7: Effective potential V n — Voo for constant curvature 
ft = 1.0. (a) On-site case, (b) Intersite case. 



For constant twist, Fig. El we find different behaviors 
in the two chain configurations. For the on-site case with 
t = 6 [Fig. |Sf a)], we see a decreasing potential well 
depth with increasing curvature at sites n = — 1 and 
n = 1, whereas the potential at n = is almost con- 
stant. This corresponds to an effective potential well for 
increasing curvature in the on-site case. In the intersite 
case [Fig. EJb)], the twist is fixed at r = 2, but in con- 
trast to the on-site case, we see that the depth of the 
potential well increases with increasing curvature at the 
center sites n = and n = 1. Therefore, bubble gen- 
eration is not found for larger curvature in the intersite 
case, corresponding to the behavior seen in Fig. [U 



C. Effective potential 

It is obvious that having initially randomly distributed 
energy along the chain, successful bubble generation 
should include, as a first stage, funneling of energy in the 
bent and twisted region. Therefore we can expect that 
only in the case when this region acts as a potential well, 
bubbling may occur. The behavior found in Figs. and 
can be qualitatively explained by the effective on-site 
potential, V n = ^2' m Jrmn which is introduced as 

n m n m 

+ J2 V " U n- 



We consider the dipole potential at given sites for both 
the on-site and the intersite case for constant curvature, 
k = 1.0, with respect to the "ground state" at n — > ±oo. 
Thus, Fig. depicts the depth of the potential well for 
constant curvature for both the on-site [Fig.CJa)] and the 
intersite [Fig. [7fb)] case . We see that in the vicinity of 
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FIG. 8: Effective potential V„ — Voo for constant twist, (a) 
On-site case, r = 6. (b) Intersite case, t — 2. 

Thus both curvature and twisting play a role in the lo- 
calized formation of precursors for denaturation bubbles 
in our model of the DNA molecule and it is clear that 
the chain configuration is important. Stronger twist in- 
creases the initiation of bubble generation in both cases 
considered. The effect of increasing curvature is different: 
In the on-site case, curvature clearly enhances bubble 
generation, but in the intersite case it slightly decreases 
the formation of bubbles in the range of k considered. 

Simulations for smaller values of the parameter J 
showed that stronger twist was required to create bub- 
bles. 



IV. ANALYTICAL APPROACH V. CONCLUSION 



Considering only the center sites n — and n = 1 in 
the intersite case, we are in effect looking at a dimer. 
Assuming that both displacements are equal, u = ui = 
u, the coupling term in the Hamiltonian, Eq. , vanishes 
and we are left with 



H = v? 



■20 



r u -i) 2 - Ju\ 



(8) 



with J = |Joi| = |Jio|- For a strong twist, r » 1, the 
last term becomes negative, due to opposite dipole ori- 
entations, corresponding to attractive interaction. The 
effective potential 



V(u) = 2 (e- u - 1) - Ju 2 

is shown as the dashed curve in Fig. [3 
Equation (|SJ may now be integrated as 



t-t = 



dw 



H -2(e- 



lf + Jw 2 



(9) 



(10) 



where u = u at the time t — i. Choosing t so large that 
w 3> 1, the power term in the square root of Eq. (|10[) . 
Jw 2 , dominates (as seen in Fig.[3Jl. Therefore, Eq. 1|10|) 
may be approximated as 



t-t 



dw 
V Jw 2 



from which we find 



u oc exp 



\fl{t - t) 



m accor- 



dance with the exponential behavior found numerically 
in Fig. 2| A similar approach can be used for the on-site 
case with identical results. 



We have shown that bubble generation in DNA-like 
models can be initiated by curvature and twisting of the 
molecular strands. 

Stronger twist facilitate bubble generation, whereas 
the effect of curvature depends on the details of the ge- 
ometry. For the on-site case, increasing curvature in- 
creases the tendency for bubble generation. Conversely, 
the intersite case decreases bubble generation for increas- 
ing curvature. 

Bubbles emerge in the region of maximal twist and 
curvature, and are found at physiological temperatures. 

Numerical results are supported by analytical approxi- 
mation. Widely used parameter values for DNA are used 
in our model. 
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